Transcriptome analysis suggested that lncRNAs regulate rapeseed seedlings in responding to drought stress by coordinating the phytohormone signal transduction pathways

The growth, yield, and seed quality of rapeseed are negatively affected by drought stress. Therefore, it is of great value to understand the molecular mechanism behind this phenomenon. In a previous study, long non-coding RNAs (lncRNAs) were found to play a key role in the response of rapeseed seedlings to drought stress. However, many questions remained unanswered. This study was the first to investigate the expression profile of lncRNAs not only under control and drought treatment, but also under the rehydration treatment. A total of 381 differentially expressed lncRNA and 10,253 differentially expressed mRNAs were identified in the comparison between drought stress and control condition. In the transition from drought stress to rehydration, 477 differentially expressed lncRNAs and 12,543 differentially expressed mRNAs were detected. After identifying the differentially expressed (DE) lncRNAs, the comprehensive lncRNAs-engaged network with the co-expressed mRNAs in leaves under control, drought and rehydration was investigated. The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of co-expressed mRNAs identified the most significant pathways related with plant hormones (expecially abscisic acid, auxin, cytokinins, and gibberellins) in the signal transduction. The genes, co-expressed with the most-enriched DE-lncRNAs, were considered as the most effective candidates in the water-loss and water-recovery processes, including protein phosphatase 2 C (PP2C), ABRE-binding factors (ABFs), and SMALL AUXIN UP-REGULATED RNAs (SAURs). In summary, these analyses clearly demonstrated that DE-lncRNAs can act as a regulatory hub in plant-water interaction by controlling phytohormone signaling pathways and provided an alternative way to explore the complex mechanisms of drought tolerance in rapeseed. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10624-4.


Introduction
Drought stress is a significant threat to global agricultural production [1].Plants respond to drought stress through morphological changes visually, such as stomata closure, leaf area reduction, and root growth promotion, as well as physiological, biochemical, and molecular mechanisms internally, including osmotic regulation, antioxidant synthesis, and upregulation of related gene expression [2].Drought stress can negatively impact plant growth, photosynthesis, respiration, and organ development, ultimately affect the crop yield and the quality of agriculture product [3].Drought resistance/tolerance mechanisms involve complex biological processes, such as gene expression, signal transduction regulation, and cellular metabolic rates with key pathways including the ABA-dependent/independent pathway and the abscisic acid signaling system [4][5][6].Under drought stress, stress signals are transmitted through various transduction pathways in plants, and some important stress signals can regulate stress-inducing genes.The expression and regulation of these related functional genes can change the morphological structure of plants or the physiological and biochemical indexes of cells to improve the drought resistance of plants [7].
In recent years, there has been a growing interest in studying the post-transcriptional level regulation to understand anti-stress mechanisms, with a focus on the function of non-coding RNAs [8,9].Long non-coding RNAs (lncRNAs) are RNA molecules longer than 200 nucleotides that do not encode proteins, which have been shown to play a role in plant resistance to drought stress [10].In Arabidopsis, the DROUGHT INDUCED lncRNA (DRIR) regulates plant response to drought stress by modulating the expression of genes involved in ABA signaling, water transport, and other stress-reducing processes [11].In rice, the lncRNA MSTRG.28732.3works with miR171 to the target genes involved in chlorophyll membrane synthesis (i.e., Os02g0662700, Os02g0663100, and Os06g0105350) in plants [12].In cotton, a large number of lncRNAs associate with ethylene, auxin, gibberellins, and cytokinins under drought stress [13].In cassava leaves and roots, 124 drought-responsive lncRNAs regulate the expression of their neighboring genes involved in hormone metabolism, transcriptional RNA regulation, and receptor kinase signaling [14].Two drought resistance-related lncRNAs were identified in tetraploid cassava, which affect the stomatal density [15].In tomato, some drought stress-related lncRNAs were found to promote the expression of target genes specifically enriched in response to the stimuli, signaling, and transporter activity [16].In Tibetan wild barley, ten lncRNAs were exclusively induced by drought stress and the lncRNA-mRNA interaction-based analysis identified the potential regulator, a serine/threonine-protein kinase SMG1 [17].However, the comprehensive survey of lncRNAs involved in drought-responsive regulation in rapeseed are still lacking; and it is unclear whether lncRNA-mRNA co-expression networks are involved in the response to drought stress.
Rapeseed (Brassica napus L., AACC, 2n = 38) is an important oilseed crop cultivated worldwide for producing edible oil, animal feed, and biodiesel.However, rapeseed is reported to be very sensitive to water deficits during the germination and seedling growth stages [18].Drought stress not only reduces the yield of rapeseed but also limits its adaptability.Therefore, drought has become one of the major limiting abiotic factors hindering the production and promotion of rapeseed, especially in China.In a previous study, we conducted a genomewide lncRNAs analysis between drought-tolerant and drought-sensitive genotypes and analyzed the possible lncRNAs' function between different genetic background under drought stress, providing some clues for understanding the drought-resistant mechanism [19].In this study, we aim to explore the lncRNA and mRNA expression divergence profiles of a drought-tolerant germplasm under control, drought plus rehydration treatments, so as to identify the key lncRNAs involved in the different water-supply environments.

Phenotype of rapeseed seedlings under different treatments
The seedlings under control (CK), drought stress (DS), and re-watering (RW) treatments exhibited different morphologies.Compared to CK, seedlings under DS showed typical dehydration symptoms, such as leaf wilting and growth retardation (Fig. 1A); however, after rehydration, seedlings appeared to recover from drought and leaves became strong and healthy again (Fig. 1A).Compared to CK, seedlings under DS had a significant reduction in fresh weight.The seedlings under RW were observed to have larger leaves.Similar with CK, the fresh weight of seedlings under RW was significantly higher than that under DS (Fig. 1B).Thus, rehydration alleviated the symptoms of drought stress in rapeseed seedlings.

Hormone contents in the leaves under different treatments
Studies have shown that phytohormones such as auxin, abscisic acid, cytokinin, gibberellins, and salicylic acid could regulate plant tolerance to drought stress [20,21].Five hormones (auxin (IAA), abscisic acid (ABA), zeatin (ZT), gibberellic acid 3 (GA3) and salicylic acid (SA)) were measured in the leaves under the CK, DS and RW treatments (Fig. 2).Compared with that under CK, the ABA content increased 37-fold under the DS treatment, the IAA content increased by 45.38%, the ZT content increased by 26.42%, the GA3 content decreased by 50.04%; moreover, the SA content decreased by 11.41%.Following the RW treatment, the ABA content nearly reached the CK level again, the IAA content was 1.25 times higher than the CK, the ZT content was restored to 1.05 times of the CK level, and the GA3 content was restored to 64.25% of the CK level, while the SA content was restored to 94.67% of the CK level.Therefore, drought stress and rehydration treatment affected the metabolism of these hormones.

RNA sequencing and data analysis
RNAs were successfully extracted from all samples (three biological replicates for each treatment) and were qualified before performing RNA sequencing (Supplementary Table 1).Clean reads were obtained by removing lowquality reads from the RNA-seq data, following standard protocols.Quality and GC-content data were then calculated from the clean data to assess the quality of the sequencing data (Supplementary Table 2).The clean datasets were then mapped to the B. napus reference genome (http://www.brassicagenome.net/databases.php).All the parameters showed the reliability of the RNA-sequencing data, and the datasets could be incredibly powerful for subsequent analysis.

Validation of sequencing data by quantitative real-time PCR (qPCR) analysis
Ten differentially expressed DE-lncRNAs were randomly selected for qRT-PCR verification.The R 2 value (R 2 > 0.9) indicated a significant correlation between the expression levels of the DE-lncRNAs quantified by FPKM and qRT-PCR results (Fig. 3).For example, the relative expression of XLOC_000799 decreased in DS vs. CK, while increased in RW vs. DS, which was consistent with the RNA-seq result (Supplementary Table 3).The real-time PCR results confirmed the expression patterns derived from transcriptome sequencing, indicating the reliability of the RNA-Seq analysis results.

Differentially expressed lncRNAs (DE-lncRNAs) and mRNAs (DE-mRNAs) under different treatment
By comparing the amount of each lncRNA and mRNA sequence from leaves in CK, DS, and RW treatments, DE-lncRNAs and DE-mRNAs were screened with a q-value threshold of < 0.05 by performing comparisons between DS vs. CK and RW vs. DS.A total of 381 DE-lncRNAs (132 down-regulated, 249 up-regulated) and 10,253 DE-mRNAs (5,377 down-regulated, 4,876 up-regulated) were identified in DS vs. CK, while 477 DE-lncRNAs (369 down-regulated, 108 up-regulated) and 12,543 DE-mRNAs (5,546 down-regulated, 6,997 up-regulated) were identified in RW vs. DS (Fig. 4).The data revealed three main findings: (1) In general, there were more DE-lncRNAs and DE-mRNAs in RW vs. DS than in DS vs. CK; (2) there were more up-regulated DE-lncRNAs in DS vs. CK than down-regulated DE-lncRNAs, while the opposite was observed in RW vs. DS; and (3) there were more down-regulated DE-mRNAs in DS vs. CK than upregulated DE-mRNAs, while the opposite was observed in RW vs. DS.

Function and pathway analysis of DE-lncRNAs based on lncRNA-mRNA co-expression network in each comparison
We used lncRNA-mRNA relationship pairs to construct an interactive network to characterize the roles and functions of DE-lncRNAs.In DS vs. CK, there were 3,493 lncRNA-mRNA pairs (including 1,423 mRNAs and 102 lncRNAs) (Supplementary Table 4).Similarly, there were   4).Gene Ontology (GO) analysis was applied on the target mRNAs to analyze their biological processes (BPs), cellular components (CCs), and molecular functions (MFs) (Fig. 5).

Deep analysis on the key DE-lncRNAs which continuously function during the water loss and water-retaining process
The Venn diagram in Fig. 7 showed the numbers of independent and overlapping DE-lncRNAs for each Fig. 5 Gene Ontology (GO) classifications of the co-expressed mRNAs of the differentially expressed lncRNAs.The mRNAs co-expressed with lncRNAs were divided into three main categories by GO analysis: biological process, molecular functions, and cellular components.The x-axis indicates the number of genes in a sub-category, and the y-axis indicates the sub-categories comparison.A total of 190 DE-lncRNAs were identified in two comparisons, of which three were up-regulated and three were down-regulated in DS vs. CK and RW vs. DS.Since dehydration and rehydration have opposite effects, the regulation direction of lncRNAs should be also opposite in these two treatments.Therefore, the remaining 184 DE-lncRNAs were the focus of subsequent analysis.Among these, 54 DE-lncRNAs were down-regulated in DS vs. CK but up-regulated in RW vs. DS, and 11 of these lncRNAs were found to have 37 partner mRNAs.Similarly, 130 DE-lncRNAs were upregulated in DS vs. CK but down-regulated in RW vs. DS, and 56 of these lncRNAs were found to have 998 partner mRNAs.These partner-mRNAs were used for subsequent Go analysis and KEGG enrichment.Based on the significance threshold of p < 0.05, these co-expressed target mRNA genes were assigned to 20 significant terms (Supplementary Table 5).In the BP category, the top five terms were oxidation-reduction process (GO:0055114), protein dephosphorylation (GO:0006470), dephosphorylation (GO:0016311), response to abiotic stimulus (GO:0009628), and embryo development (GO:0009790).In the MF category, the top five terms were nucleic acid binding transcription factor activity (GO:0001071),  6 KEGG pathway analysis.Top 20 pathways for the co-expressed mRNAs of the differentially expressed lncRNAs.The y-axis corresponds to the KEGG pathway with a q-value ≤ 0.05, and the x-axis shows the enrichment ratio between the number of differentially expressed genes and all UniGenes enriched in a particular pathway.The color of the dot represents q value, and the size of the dot represents the number of differentially expressed genes mapped to the reference pathways.(A) KEGG pathway classification of the mRNAs co-expressed with DE-lncRNAs in DS vs. CK.(B) KEGG pathway classification of the mRNAs co-expressed with DE-lncRNAs in RW vs. DS sequence-specific DNA binding transcription factor activity (GO:0003700), sequence-specific DNA binding (GO:0043565), phosphoprotein phosphatase activity (GO:0004721), and protein serine/threonine phosphatase activity (GO:0004722).Notably, CCAAT-binding factor complex (GO:0016602) was the only enriched term in the CCs category.
KEGG pathway analysis was separately conducted in two comparison groups to understand the biological pathways in water-deficient or water-sufficient environments (Supplementary Table 6).Three pathways, including plant hormone signal transduction (35 genes in DS vs. CK and 36 genes in RW vs. DS), carbon metabolism (27 genes in DS vs. CK and 30 genes in RW vs. DS), and glycolysis/gluconeogenesis (17 genes in DS vs. CK and 17 genes in RW vs. DS), were significantly enriched in both comparison groups with the largest number of genes.
Because too many GO analysis results were enriched in the comparison, we used the KEGG results to focus on the major melatonin-related pathways and pivotal genes, which would emphasize the main melatonin function.
Due to the excessive GO analysis results in the two comparison series (DS vs. CK and RW vs. DS), we prioritized the most important drought-related pathways and key genes based on the KEGG results.Interestingly, KEGG analysis revealed that most genes in the plant hormone signal transduction pathway were enriched.To further investigate the mechanism of lncRNAs functioning under drought stress, we focused on analyzing the lncRNAs co-expressed with genes involved in plant signal transduction pathways.

Analysis of DE-lncRNAs involved in plant hormone signal transduction pathways
Since over 15% of co-expressed mRNAs of DE-lncRNAs were enriched in the pathway of plant hormone signal transduction in both comparisons, this pathway was selected as the focus.In DS vs. CK, the co-expression network under this pathway contained 108 matched lncRNA-mRNA pairs, including 24 lncRNAs and 35 mRNAs (Fig. 8A and Supplementary Table 7).The target genes of 24 DE-lncRNAs consisted of 1 down-regulated mRNA and 34 up-regulated mRNAs.In RW vs. DS, the co-expression network of plant hormone signal transduction contained 157 matched lncRNA-mRNA pairs, including 41 lncRNAs and 36 mRNAs (Fig. 8B and Supplementary Table 7).The target genes of the 41 DE-lncRNAs consisted of 35 down-regulated mRNAs and 1 up-regulated mRNA.As shown in Fig. 9, the expression trends or levels of these lncRNAs co-expressed with auxin, cytokinins, gibberellins, abscisic acid, ethylene, and salicylic acid signal transduction genes were completely different under drought and rewater treatment.Within this pathway, the target genes involved in abscisic acid (ABA) signal transduction were dominant (65.7% and 66.7% in the two comparisons).The genes with the highest co-expression frequency included type 2 C protein phosphatases (PP2Cs) and abscisic acidresponsive transcription factors (ABFs) in the ABA metabolism pathway, as well as small auxin upregulated RNAs (SAURs) in the IAA metabolism pathway.The expression levels of lncRNAs co-expressed with some genes involved in plant hormone signal transduction (PP2Cs, ABFs, and SAURs) were measured using qRT-PCR, and the results of qRT-PCR were consistent with the RNA-seq data, indicating that these lncRNAs were key genes in response to drought and rehydration (Fig. 10).

Discussion
Abiotic stress seriously affects agricultural development, and lncRNAs have been shown to take important parts in regulate crops' response to the environments at the molecular level [22].Water supply is a vital environmental factor in agriculture, and excessive water or insufficient water (drought) supply can both have a negative impact on crops [23].Studies have shown that lncRNAs function under drought stress, but their molecular mechanism is not yet well understood.In previous studies, we only compared rehydration treatment with drought stress and identified lncRNAs related to rehydration as involved in plant hormone signal transduction.So, what pathways and functions do lncRNAs induced by drought stress participate in?In this study, libraries of both lncRNA and mRNA were constructed under normal conditions (CK), drought stress (DS), and re-watering (RW).The contents of five phytohormones were measured under different treatment conditions.Our data supported the comparisons of the whole process, from normal to drought and then to rehydration, at the physiological level and the molecular level.We presented a comprehensive analysis to uncover the crucial role of important lncRNAs in the process of water loss and rehydration.

Phenotypic changes of rapeseed seedlings during waterlosing and rehydration
The phenotypic and physical parameters indicated significant differences in the growth status of rapeseed seedlings under various water supply conditions.The DS treatment decreased the seedling weight significantly and the RW treatment recovered it to a similar level as CK, indicating that the temporary damage from drought is reversible and repairable in rapeseed seedlings, which is according with Xu et al. [24].However, the seedling phenotype uder RW was not consistant with that under DS, implying the molecular-level activities happened during the water loss and water supply enviroments.

Genes related with the phytohormones signal transduction during drought stress
Although 190 co-expressing DE-lncRNAs were identified as being related to many biological processes, the number of lncRNAs involved in plant hormone signal transduction was the largest.In this pathway, target genes of DE-lncRNAs were mostly involved in abscisic acid, auxin, cytokinin, and gibberellins signaling pathways in both comparison groups.
Abscisic acid (ABA) is a key phytohormone that is essential in regulating various growth and metabolic processes in response to biotic and abiotic stresses, especially drought stress [25].The increase in ABA concentration in the soil solution under drought stress stimulates ABA signaling in roots and stems, regulating the water status of plants [26,27].The ABA signaling system includes the ABA receptor proteins PYR/PYL/RCAR, the positive regulator SnRK2 (SNF1 associated protein kinase 2), and the negative regulator 2 C protein phosphatases (PP2C), as well as their downstream targets.These three components-PYR/PYL/RCAR, SnRK2, and PP2C-combine as a Fig. 9 Expression levels of transcripts involved in phytohormone signal transduction pathways in DS vs. CK and RW vs. DS dual negative regulatory system to regulate ABA signaling and its downstream responses [28,29].
However, some studies have shown that subfamily A PP2Cs in Arabidopsis and other plants negatively regulates ABA and stress signaling pathways [30][31][32].On the other hand, BdPP2CA6, a subfamily A PP2C in Brachypodium distachyon, was found to be a positive regulator of ABA and stress signaling pathways [33].Similarly, the expression of SiPP2C10 in foxtail millet has been upregulated under drought stress conditions [34].This finding was consistent with the results of our study, leading us to speculate that these lncRNAs co-expressed with PP2C may be closely related to drought resistance genes.
In this study, we found that both the two core components of the ABA signal transduction pathway and the downstream central ABF component were regulated by drought stress.When rapeseed seedlings were exposed to drought stress during early growth stages, ABA content was significantly increased.The expression of lncRNAs co-expressed with PP2C, SnRK2 and ABF coding genes was upregulated by high concentration of ABA, resulting in enhanced ABA signal, which may increase the inhibition of ABA on the growth of rapeseed seedlings (Supplementary Figure S1A).After rehydration, ABA content almost recovered to the control level, and the expression of lncRNAs co-expressed with PP2C, SnRK2, and ABF was down-regulated, resulting in the weakening of the ABA signal, which alleviated the inhibitory effect of ABA on the growth of rapeseed seedlings (Supplementary Figure S1B).The analysis of the lncRNA-mRNA co-expression network showed that the genes encoding PP2C, SnRK2, and ABF had similar expression patterns during both the drought and rewatering processes.However, the lncRNAs co-expressed with these genes had not only different expression patterns, but also different numbers, revealing that rehydration is not just a simple process of Fig. 10 Confirmation of the expression patterns of key lncRNAs under drought and rehydration conditions using quantitative RT-PCR restoring drought stress, but also involves lncRNAs coexpression with mRNAs in the regulation of a diverse range of recovery processes.
Auxin (IAA) is a major regulatory signal for plant growth and development and negatively regulates plant drought resistance [35,36].The IAA early response genes, such as GH3, Aux/IAA, and SAUR, play a role in the IAA signaling pathway [37,38].The Aux/IAA protein functions as a transcriptional inhibitor in the IAA signal transduction pathway [39].GH3 encodes auxinconjugating enzymes that are involved in stress responses by controlling the levels of active auxin through negative feedback [40].The SAUR39 gene of rice acts as a negative regulator of auxin synthesis and transport [41].The co-expression network of lncRNA-mRNA, Aux/IAA, GH3, and SAUR genes in the auxin signal transduction pathway was analyzed in response to drought stress and rehydration.
IAA content increased significantly after drought stress, the signaling pathway was activated, and lncRNAs co-expressed with Aux/IAA, GH3, and SAUR were upregulated, indicating that stress restricted the growth of rape seedlings by inhibiting IAA signal transduction (Supplementary Figure S2A).After rehydration, IAA content decreased, resulting in decreased expression of lncRNAs co-expressed with Aux/IAA, GH3, and SAUR genes, which accelerated the vegetative growth by stimulating tissue elongation and cell expansion (Supplementary Figure S2B).
Cytokinins (CTKs) not only manage plant growth and development but also play a role in mediating plant tolerance to drought stress [42].The accumulation of CTKs can have both positive and negative effects on plant drought resistance [43].Arabidopsis seedlings use a twocomponent signaling system (TCS) to regulate cytokinin signaling, which consists of sensor histidine kinases (AHKs), histidine phosphate transfer proteins (AHPs), and response regulators (ARRs) [44].Cytokinin Response 1 (CRE1) was found to negatively regulate osmotic pressure in the presence of CTKs [45].In Arabidopsis, AHP2, AHP3, and AHP5 serve as redundant negative regulators in response to drought stress [46].The expression B-type cytokinin response regulators ARR1, ARR10, and ARR12 in Arabidopsis were inhibited by drought stress, indicating that they negatively regulate plant drought tolerance [47].
The content of zeatin, the main natural active component of CTKs, increased under drought stress, which caused the up-regulation of lncRNAs co-expressed with CRE1, AHP and B-ARR.Compared with drought treatment, zeatin content decreased after rehydration treatment, which led to downregulation of lncRNAs coexpressed with CRE1, AHP and B-ARR (Supplementary Figure S3).Results of the study indicated that drought stress inhibited cytokinin signal transduction, leading to reduced growth of rapeseed seedlings.Conversely, after rehydration, an increase in cytokinin signal transduction resulted in improved growth of the seedlings.Additionally, CRE1, AHP, and B-ARR were found to play a role as negative regulators in the growth of rapeseed seedlings under drought stress.
Interestingly, the co-expression of lncRNA-mRNA revealed that while CRE1 (BnaC07g11340D), AHP (BnaC01g31940D), and B-ARR (BnaA01g17750D) were involved in cytokinin signal transduction during both drought and rehydration, and the lncRNAs co-expressed with these genes in the two treatments differed.It suggested that rehydration was not merely a simple drought recovery process, but that lncRNAs also played a role in this process.
Gibberellins (GAs) not only regulate seed germination, stem elongation, and flower development but also participate in the regulation of abiotic stress [48,49].There are many kinds of gibberellins, among which GA1, GA3, GA4, and GA7 have the highest biological activity [49].The synthesis of gibberellin is inhibited by drought stress [50].Research has shown that reducing GAs levels can improve drought resistance and help plants overcome water deficits through drought avoidance [50,51].The sensing of gibberellins (GAs) is mediated by GID1 (GA-INSENSITIVE DWARF1), a receptor like a hormone-sensitive lipase [52].When GA binds to GID1s, it stimulates interaction between GID1s and the growthinhibiting DELLA [53].
In our experiment, the content of GA3 in leaves decreased significantly under drought stress.After re-watering treatment, the content of GA3 increased significantly compared with that after drought treatment, although the level of GA3 was not up to the level of control.During drought stress, the expression of the lncRNA (XLOC_042894) that co-expressed with BnaCnng55170D (encoding GID1) was up-regulated (Supplementary Figure S4A).After rehydration, the expression of lncRNAs (XLOC_098397, XLOC_038342, XLOC_015081, XLOC_071559, XLOC_100682, and XLOC_042894) that co-expressed with BnaA07g19530D and BnaCnng55170D (encoding GID1) were down-regulated (Supplementary Figure S4B).The results showed that drought stress reduced the content of GA3 in rapeseed seedlings, and GID1 did not bind with low concentrations of GA3, which made the DELLA protein bind to the gibberellin response gene and inhibit its activity, thus inhibiting the growth of rapeseed seedlings.Rehydration inhibited the repression of DELLA on GA signaling, resulting in improved growth of rapeseed seedlings.
In conclusion, the number of lncRNAs co expressed with the genes encoding PP2C and ABF is the highest in the ABA signaling pathway.Moreover, compared to drought stress, rehydration treatment mobilized more lncRNAs to participate in ABA, IAA, CTKs, and GAs signal transduction, helping rapeseed seedlings recover growth.

The proposal model showing LncRNAs' function in response to drought stress and rehydration in rapeseed seedings
By focusing on DE-lncRNAs, we proposed possible mechanisms.Under drought stress, certain specific lncRNAs co-expressed with genes were responsible for the corresponding hormone's signal transduction, regulating the expression of downstream genes and affecting seedlings growth.For instance, under drought stress, 21, 2, and 15 lncRNAs were up-regulated and co-expressed with PP2C, SnRK2, and ABF, respectively, promoting ABA signal transduction and intensifying the inhibitory effect of ABA on the rapeseed seedlings.After rehydration, 35, 5, and 24 lncRNAs co-expressed with PP2C, SnRK2, and ABF were down-regulated, releasing the inhibition, and allowing the seedlings to continue growing.It should be noted that ABA, IAA, and GAs are the major hormones involved in directing seedling growth (Fig. 11).

Plant materials, growth conditions, and treatments
A drought-tolerant B. napus germplasm Q2, obtained from the Oil Crops Research Institute (OCRI) of the Chinese Academy of Agricultural Sciences (CAAS), was used in this study.The high resistance of Q2 to drought was confirmed in a previous study [54].This experiment was conducted in a controlled greenhouse (20 °C, 16 h/8 h L/d, 65% humidity).Q2 seeds were sterilized with 0.1% HgCl 2 for 10 min, and rinsed three times in the sterile water for 20 min.Seeds were put on the moist filter paper in the petri dishes for germinating at 25 °C.After 7 days, uniform seedlings were randomly selected and transplanted into plastic pots (6.5 cm × 4.5 cm), filled with a mixture of 70 g of soil, vermiculite, and sand (soil: vermiculite: sand = 2:1:1, v/v/v).The vermiculite mixture had a field capacity (FC) of 45.21% determined by the methods of Wilcox [55] and Duan et al. [56].The plants were watered normally from germination until the three-leaf stage (at 18 days), followed by different treatments: (1) normal watering for 8 days and maintaining 75% FC (set as control, CK); (2) 8-day drought treatment (no watering), decreasing the FC to 35% (set as drought stress, DS) [57]; and (3) 7-day drought treatment followed by 1-day of rehydration (adding water), causing the FC to recover to 75% (set as re-watering, RW).This experiment sets as a completely randomized design with 3 replicates.After each treatment, the third leaf (from the top of the seedlings) was cut into five parts and mixed to form one sample for each replicate of each treatment.All samples (9 in total) were flash frozen in liquid nitrogen and stored at -80 °C prior to RNA extraction.

Determination of physiological parameters and hormone contents
The uniform-growing seedlings were selected under the CK, DS, and RW treatments.After recording the fresh weight of the plants, seedlings were divided into shoots and roots.The content of five endogenous plant hormones (IAA, ABA, ZT, GA3, and SA) in the leaves of rapeseed seedlings under different treatments was measured via HPLC-MS/MS analysis [58,59].

RNA extraction and testing
The total RNA was extracted from each sample according to the manufacturer's instruction using Trizol reagent (Invitrogen, Burlington, ON, Canada).RNA quality was assessed using several methods: firstly, the RNA degradation and contamination were monitored using the 1% agarose gel; secondly, RNA purity was detected by the Nanophotometer spectrophotometer (IMPLEN, CA, USA).Thirdly, RNA concentration was measured using the Qubit RNA Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, CA, USA).Finally, RNA integrity was assessed using the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, USA).

Library construction, sequencing, and mapping to the reference genome
The library preparation and deep sequencing were performed by the Novogene Bioinformatics Technology Cooperation (Beijing, China).A total amount of 3 µg RNA per sample was used as input material for the RNA sample preparations.Firstly, ribosomal RNA was removed by Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, USA), and rRNA free residue was cleaned up by ethanol precipitation.Subsequently, sequencing libraries were generated using the rRNA-depleted RNA by NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following manufacturer's recommendations.Briefly, fragmentation was carried out using divalent cations under elevated temperature in NEBNext First Strand Synthesis Reaction Buffer (5X).First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH-).Second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H.In the reaction buffer, dNTPs with dTTP were replaced by dUTP.Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities.After adenylation of 3' ends of DNA fragments, NEB-Next Adaptor with hairpin loop structure were ligated to prepare for hybridization.In order to select cDNA fragments of preferentially 150 ~ 200 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA).Then 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR.Then PCR was performed with Phusion High-Fidelity DNA polymerase, Universal PCR primers and Index (X) Primer.At last, products were purified (AMPure XP system) and library quality was assessed on the Agilent Bioanalyzer 2100 system.Each sample was processed individually, and named as CKQ2-1, CKQ2-2, CKQ2-3, DSQ2-1, DSQ2-2, DSQ2-3, RWQ2-1, RWQ2-2, and RWQ2-3, resulting in a total of 9 libraries.All libraries were sequenced on the Illumina Hiseq 2500 platform, generating 125 bp paired-end reads.The original data of this experiment has been uploaded to the NCBI database, and the relevant accession numbers is PRJNA876031.Quality parameters such as Q20, Q30, and GC content of the clean data were calculated.Clean reads from each sample library were mapped to the reference database of Brassica genomes (https://appliedbioinformatics.com.au/gb2/gbrowse/BnapusPan/) using Bowtie v2.0.6.

Identification of lncRNAs
Before screening, the Cuffmerge software was used to merge transcripts spliced from each sample and remove transcripts with an uncertain chain direction to obtain complete transcriptome information for this sequencing.Then, we screened the combined transcript set for lncRNA using the following steps: (1) Select transcripts with a number of exons ≥ 2; (2) Select transcripts with a length > 200 bp; (3) Screen out transcripts that overlap the exon region of the database annotation using Cuffcompare software and include lncRNA that overlap the exon region of the transcript splicing in the database as database annotation lncRNA for subsequent analysis; (4) Calculate the expression amount of each transcript using Cuffquant, and select transcripts with FPKM ≥ 0.5; (5) For spliced transcripts, we employed bioinformatics tools such as CPAT, CNCI, PfamScan, and phyloCSF to comprehensively assess the coding potential of transcripts for screening, taking the intersection of predicted transcripts without coding potential in the analyzed software results.Subsequently, we identified transcripts that were predicted to have coding potential by at least one of the coding potential prediction software as TUCP (transcripts of uncertain coding potential) for this analysis.Among these transcripts, there may be a subset of lncRNA with certain coding potential, and therefore, we included them as a separate category of transcripts for subsequent analysis [60].

Quantification of gene expression level, and gene expression profiles
Firstly, the Cuffdiff software was used to calculate the FPKMs (fragments per kilo-base of exon per million fragments mapped) of both lncRNAs and mRNAs in each sample.Then, the FPKM method was used to calculate the expression level of each transcript.Finally, the average value of the three replicates from the same treatment was recorded as the expression intensity of specific genes [61,62].Differentially expressed lncRNAs and mRNAs were identified using the Cufflinks software with a less stringent threshold of q-value ≤ 0.05 for the significant gene expression comparison between samples.The q-value, which is the corrected p-value, was used to estimate the false discovery rate (FDR) in multiple comparisons [63].The expressions of lncRNAs and mRNAs were then compared under different treatments (DS/CK, RW/ DS).

Construction of the lncRNA-mRNA co-expression network
Since the function of lncRNAs is not well-defined, a commonly used method to predict their functional mechanism is based on co-expressed mRNAs [64].Moreover, previous studies have shown that some lncRNAs can form "lncRNA-mRNA pairs" with nearby protein-coding genes, which can affect their function [65].To explore the co-expression relationship between lncRNAs and mRNAs, we constructed a co-expression network using the method described by Wang et al. [66].Pearson's correlation coefficient (PCC) and p-value were calculated for each lncRNA-mRNA pair, and we selected pairs with |PCC value| ≥ 0.95 and p < 0.05 to construct the network.The network was visualized using Cytoscape (v3.7.1; https://cytoscape.org/).

GO and KEGG enrichment analysis
We performed GO enrichment analysis on the target genes that were co-expressed with differentially expressed lncRNAs using the GOseq R package [67].GO terms were used to classify the genes into different categories, such as biological processes, molecular functions, and cellular components.The threshold of p-value < 0.05 was used to determine significantly enriched GO terms.
The KOBAS software was used to detect the statistical enrichment of differentially expressed lncRNAs target genes in the KEGG pathway (http://www.genome.ad.jp/ kegg/) [68].We considered a p-value < 0.05 to indicate statistically significant differences between samples.

Validation of DE-lncRNAs through real-time quantitative PCR
To obtain the first-strand cDNA, the total RNA of each sample was treated with RNase-free DNase, and then the RevertAid First-Strand cDNA Synthesis Kit (Fermentas, USA) was used.For real-time PCR, a 10 µl reaction system was prepared on the ABI 700 Real-time PCR platform containing approximately 0.5 ng of cDNA, 2.5 µl of 1.2 µM mixture of forward and reverse primers, and 5 µl of master mix, as directed by the SYBR Green PCR Master Mix system (Takara Co. Ltd., Japan).The PCR amplification conditions were set as follows: one cycle of 95 °C for 30 s, followed by 40 cycles of 95 °C for 5 s, and 60 °C for 30 s.We used Primer Premier 5 software to design and analyze PCR primers for the detection of candidate lncRNAs, which were listed in Supplementary Table 8.In our experiment, random primers were used for reverse transcription.For each experiment, three independent biological replicates were performed, and for each sample in the RT-qPCR reaction, three biological and three technical replicates were performed.

Conclusion
In this study, seedlings of a drought-tolerant rapeseed germplasm were subjected to drought and rehydration treatments.The transcriptomes of lncRNAs and mRNAs were analyzed through sequencing of the leaves under CK, DS, and RW conditions.A total of 184 DE-lncRNAs were found to be consistently expressed in the two comparison groups of DS vs. CK and RW vs. DS.A lncRNA-mRNA network was established to understand the role of lncRNAs in responding to drought stress and rehydration.The results showed that the plant hormone signal transduction pathway was the most significantly enriched in both DS vs. CK and RW vs. DS according to the transcriptome assay.Several enriched candidate mRNAs affecting phytohormone function were identified, and they played a role in drought stress tolerance in rapeseed seedlings through interaction with related lncRNAs.The study revealed a multiple plant hormone signaling mediated network regulating drought resistance in rapeseed through the lncRNA-mRNA network, with lncRNAs induced by plant hormone signals involved in ABA, IAA, CTKs, and GAs signal transduction pathways.This is the first discovery of the lncRNA-mRNA network involved in rapeseed under drought stress and drought-rehydration process.The findings provided new insights into lncRNAs in response to water loss & supply, enriched the related plant hormone signal transduction pathways and key genes for understanding drought tolerance mechanisms in plants.Manipulation of the candidate lncRNAs or the co-expressing mRNAs may enhance crop drought tolerance in the future.

Fig. 1 Fig. 2 Fig. 4
Fig. 1 Phenotype of seedlings under different treatments.(A) Phenotypes of seedlings under CK, DS, and RW treatments.Bar = 1 cm.(B) Comparisons of plant height among the treatments.(C) Comparisons of fresh weight among the treatments.The experiments were repeated three times and vertical bars indicated standard errors.CK = control; DS = drought stress; RW = re-watering.The same letter indicates no significant difference and the different letters indicate a significant difference at a significance level of 5%

Fig. 3
Fig. 3 Validation of the expression levels of lncRNAs using the real-time quantitative polymerase chain reaction (RT-qPCR).The x-axis indicates the log 2 (Fold change) as measured by RT-qPCR.The y-axis indicates the log 2 (Fold change) as measured by RNA sequencing (RNA-seq).The squared of the Pearson's correlation coefficient of relative expression measured by RNA-seq and RT-qPCR was 0.91152

Fig. 7
Fig. 7 Venn diagram showing the number of unique and common DE-lncRNAs in two comparison groups

Fig. 8
Fig. 8 LncRNA-mRNA-network analysis of the plant hormone signal transduction.The circle and rectangle nodes represent lncRNAs and protein-coding genes, respectively.The up-regulated and down-regulated nodes are separately colored in red and green.Edges show regulatory interactions among nodes.(A) 24 DE-lncRNAs interacted with 35 mRNAs in the meaningful "plant hormone signal transduction" in DS vs. CK.(B) 41 DE-lncRNAs interacted with 36 mRNAs in the meaningful "plant hormone signal transduction" in RW vs. DS

Fig. 11
Fig.11The model for the function of lncRNAs in response to drought stress in rapeseed seedlings.The up-pointing red arrows mean that the candidate lncRNAs are up-regulated; the down-pointing green arrows mean that the candidate lncRNAs are down-regulated.The main affected pathways were shown in the boxes.The number of important lncRNAs affecting the corresponding pathway is indicated in brackets.ABA = abscisic acid; IAA = auxin; CTKs = cytokinins; GAs = gibberellins; ZT = zeatin; GA3 = gibberellic acid 3